!> Background grid list (2-dimensional)
!> 背景网格链表法 (2 维)
module linked_list_nnps_m

    use list_type_m, only: list_t, list_init => list_t_init, list_push => list_t_push, merge => list_t_merge
    use twin_list_type_m, only: twin_list_t, twin_init => twin_list_t_init, twin_push => twin_list_t_push
    implicit none
    private

    public :: query, build_list_from_loc
    real, parameter :: sqrt_eps = sqrt(epsilon(1.0))

contains

    !> Build background grid (build storage structure)
    !> 建立背景网格 (构建存储结构)
    pure subroutine setup_list(hsml, min, max, estimated_capacity, list)
        real, intent(in) :: hsml                                !! smoothing length
                                                                !! 光滑长度
        real, intent(in) :: min(2)                              !! minimum value
                                                                !! 最小值
        real, intent(in) :: max(2)                              !! maximum value
                                                                !! 最大值
        integer, intent(in) :: estimated_capacity               !! estimated capacity
                                                                !! 预估表存储元素量
        type(list_t), intent(out), allocatable :: list(:, :)    !! grid
                                                                !! 网格
        integer :: i, j

        i = int((max(1) - min(1))/hsml) + 1
        j = int((max(2) - min(2))/hsml) + 1

        allocate (list(i, j))
        call list_init(list, estimated_capacity)

    end subroutine setup_list

    !> Deploy particle to background grid (mapping)
    !> 部署粒子至背景网格 (映射粒子，mapping)
    pure subroutine mapping_to_list(i, loc, min, hsml, list)
        integer, intent(in) :: i                    !! particle index
                                                    !! 粒子编号
        real, intent(in) :: loc(2)                  !! particle location
                                                    !! 粒子坐标
        real, intent(in) :: min(2)                  !! minimum value
                                                    !! 最小值
        real, intent(in) :: hsml                    !! smoothing length
                                                    !! 光滑长度
        type(list_t), intent(inout) :: list(:, :)   !! grid
                                                    !! 背景网格
        integer :: index(2)

        index(:) = int((loc(:) - min(:))/hsml) + 1

        call list_push(list(index(1), index(2)), i)

    end subroutine mapping_to_list

    !> Construct background grid from particle location
    !> 从粒子位置构建背景网格
    pure subroutine build_list_from_loc(loc, hsml, estimated_capacity, list)
        real, intent(in) :: loc(:, :)                           !! particle location
                                                                !! 粒子坐标
        real, intent(in) :: hsml                                !! smoothing length
                                                                !! 光滑长度
        integer, intent(in) :: estimated_capacity               !! estimated capacity
                                                                !! 预估表存储元素量
        type(list_t), intent(out), allocatable :: list(:, :)    !! grid
                                                                !! 背景网格
        integer :: i

        associate (min => minval(loc, dim=2) - sqrt_eps)
            call setup_list(hsml, min, maxval(loc, dim=2), estimated_capacity, list)
            do concurrent(i=1:size(loc, 2))
                call mapping_to_list(i, loc(:, i), min, hsml, list)
            end do
        end associate

    end subroutine build_list_from_loc

    !> Query twin particle list
    !> 搜索 (query)
    !> @note L 型查找方式，借鉴 SPHysics 的方案
    pure subroutine query(loc, list, hsml, estimated_capacity, twin_list)
        real, intent(in) :: loc(:, :)                   !! particle location
                                                        !! 粒子坐标
        type(list_t), intent(in) :: list(:, :)          !! grid
                                                        !! 网格
        real, intent(in) :: hsml                        !! smoothing length
                                                        !! 光滑长度
        integer, intent(in) :: estimated_capacity       !! estimated capacity
                                                        !! 预计粒子对数量
        type(twin_list_t), intent(inout) :: twin_list   !! twin list
                                                        !! 粒子对
        type(list_t) :: found
        integer :: i, j, ik, m, n, tmp(2)

        call twin_init(twin_list, estimated_capacity, .true.)
        m = size(list, 1); n = size(list, 2)
        ! 仅处理中间区域 (\[i \in [2, size(list, 1)-1],  j \in [2, size(list, 2)-1]\])
        do concurrent(j=2:n - 1, i=2:m - 1)
            call list_init(found, 64)
            if (list(i, j)%log == 0) cycle
            ! 查询项链的 L 型网格
            call merge(found, list(i - 1, j - 1))
            call merge(found, list(i, j - 1))
            call merge(found, list(i + 1, j - 1))
            call merge(found, list(i + 1, j))

            call find_nearby_particles(list(i, j), found, hsml, twin_list)
        end do

        ! i = [1, m], j = 2, n
        tmp = [1, m]
        do ik = 1, 2
            i = tmp(ik)
            do j = 2, n
                call list_init(found, 64)
                if (list(i, j)%log == 0) cycle
                call merge(found, list(i, j - 1))

                call find_nearby_particles(list(i, j), found, hsml, twin_list)
            end do
        end do

        ! i = 2, m - 1, j = [1, n]
        tmp = [1, n]
        do ik = 1, 2
            j = tmp(ik)
            do i = 2, m - 1
                call list_init(found, 64)
                if (list(i, j)%log == 0) cycle
                call merge(found, list(i + 1, j))

                call find_nearby_particles(list(i, j), found, hsml, twin_list)
            end do
        end do

    contains

        !> Find nearby particles
        !> 查找邻近粒子 (find_nearby_particles)
        !> @note 简化代码
        pure subroutine find_nearby_particles(list, found, hsml, twin_list)
            type(list_t), intent(in) :: list                !! grid
                                                            !! 表
            type(list_t), intent(in) :: found               !! nearby particles
                                                            !! 表
            real, intent(in) :: hsml                        !! smoothing length
                                                            !! 光滑长度
            type(twin_list_t), intent(inout) :: twin_list   !! twin list
                                                            !! 粒子对
            integer :: ii, ij

            do ii = 1, list%log
                do ij = 1, list%log
                    if (ii == ij) cycle
                    if (norm2(loc(:, list%items(ii)) - loc(:, list%items(ij))) < hsml) then
                        call twin_push(twin_list, [list%items(ii), list%items(ij)])
                    end if
                end do

                do ij = 1, found%log
                    if (norm2(loc(:, list%items(ii)) - loc(:, found%items(ij))) < hsml) then
                        call twin_push(twin_list, [list%items(ii), found%items(ij)])
                    end if
                end do
            end do

        end subroutine find_nearby_particles

    end subroutine query

end module linked_list_nnps_m
